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Abstract  In  a  previous  paper  [1],  an  approach  to  celestial  navigation  was  presented  that 
allows  a  vessel’s  latitude,  longitude,  course,  and  speed  to  be  simultaneously  estimated.  The 
development  assumed  that  the  vessel  was  sailing  a  rhumb-line  track  at  a  constant  speed. 
In  this  paper,  the  approach  is  extended  to  cover  the  case  where  observations  are  taken 
over  several  legs  of  a  voyage.  A  single  least-squares  solution  for  all  sailing  parameters  is 
developed.  The  algorithm  presented  in  this  paper  is  not  limited  to  celestial  navigation,  but 
is  applicable  to  any  situation  in  which  observations  over  several  legs  of  a  voyage  are  to  be 
combined  in  a  single  solution. 

Key  Words  celestial  navigation,  celestial  fix,  motion  of  observer 


Introduction 

A  development  of  celestial  navigation  was  presented  in  [1]  in  which  the  motion  of  the  observer  was 
included  as  an  essential  part  of  the  basic  mathematics.  This  allows  celestial  observations  to  be  used 
for  the  determination  of  not  just  latitude  and  longitude  at  a  specific  instant  (the  traditional  fix),  but 
also  course  and  speed — providing,  of  course,  that  a  sufficient  number  of  observations  are  available, 
suitably  distributed  in  azimuth  and  time.  The  development  in  [1]  assumes  that  the  observer’s  vessel 
is  sailing  a  rhumb-line  track  at  constant  speed  during  the  observations.  For  altitude  sights  made  with 
a  hand-held  sextant  (typical  accuracy  ±1  arcmin),  at  least  eight  observations,  spread  over  several 
hours,  are  generally  needed  for  reliable  course  and  speed  estimates.  However,  the  time  span  over 
which  the  observations  need  to  be  made  would  be  less  if  their  accuracy  were  better.  For  example,  an 
automated  shipboard  star  tracker  might  make  the  necessary  observations  and  obtain  a  good  solution 
in  a  matter  of  minutes. 

The  algorithm  described  in  [1]  has  been  implemented  in  software  specifically  designed  for  U.S. 
Navy  operational  use.  Navy  ships  are  subject  to  changes  in  orders  en  route,  and  are  often  involved 
in  maneuvers.  The  software  was  therefore  required  to  handle  an  arbitrary  number  of  changes  to 
course  and  speed.  Even  on  an  uncomplicated  port-to-port  voyage,  the  great-circle  route  would 
usually  be  approximated  by  a  series  of  rhumb  lines.  The  possibility  that  a  celestial  fix  might  involve 
observations  that  span  two  or  more  voyage  legs  had  to  be  considered. 

In  chart-based  navigation,  lines  of  position  (whether  celestial  or  terrestrial)  can  be  easily  ad¬ 
vanced  or  retired  over  changes  to  course  and  speed  (see  Bowditch  [2],  pp.  130-132).  The  assumption 
is  simply  made  that  the  difference  between  a  vessel’s  true  position  and  its  estimated  position — the 
error  in  position — remains  constant,  in  azimuth  and  distance,  as  the  vessel  moves,  despite  course 
and  speed  changes.  The  procedure  for  advancing  a  line  of  position  is  thus  no  more  complicated  when 
several  voyage  legs  are  involved  than  when  there  is  only  a  single  leg  to  consider.  A  popular  mathe¬ 
matical  approach  to  celestial  sight  reduction  [3,4]  is  a  direct  mathematical  translation  of  chart-based 
navigation,  and  it  adopts  the  same  assumption  for  combining  observations  spanning  multiple  legs. 

However,  as  pointed  out  in  [5],  even  when  only  one  leg  is  involved,  this  assumption  cannot  be 
rigorously  true,  and  it  is  often  not  even  approximately  true.  A  vessel’s  estimated  course  and  speed 
over  bottom  may  be  considerably  in  error  due  to  components  of  current  and  wind  not  taken  into 
account.  In  practice,  then,  the  error  in  position  changes  with  time.  The  situation  is  complicated 
further  by  changes  to  direction  or  speed,  which  modify  the  net  force  of  the  wind  on  the  vessel.  So 
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not  only  is  the  basic  assumption  (constant  positional  error)  dubious,  but  so  is  its  first  derivative 
(constant  rate  of  change  of  positional  error)  when  multiple  legs  are  involved. 

The  development  presented  in  [1]  assumed  at  the  outset  that  the  observer’s  estimated  course  and 
speed  are  probably  not  accurate,  so  that  the  error  in  position  is  not  constant.  That  algorithm  allows 
for  a  simultaneous  solution  for  latitude,  longitude,  course,  and  speed.  However,  the  algorithm 
is  based  on  the  geometry  of  a  single  rhumb-line  leg — that  is,  the  course  and  speed  are  assumed 
constant.  This  paper  explores  how  the  approach  in  [1]  can  be  generalized  for  the  multiple-leg  case. 
The  multiple-leg  algorithm  presented  here,  although  developed  for  rhumb-line  tracks  and  celestial 
navigation,  is  applicable  to  any  kind  of  track  and  any  set  of  observations  in  which  a  single  observable 
is  a  known  function  of  the  observer’s  instantaneous  latitude  and  longitude. 


Single-Leg  Solution 

The  algorithm  presented  in  [1]  is  a  least-squares  solution  for  a  vessel’s  latitude  and  longitude 
at  a  specific  time,  along  with  its  course  and  speed  (assumed  constant).  We  assume  that  initially  we 
know  the  vessel’s  track  reasonably  well.  Errors  of  several  degrees  in  position  or  course  or  several 
knots  in  speed  are  routinely  handled,  and  larger  errors  usually  only  require  additional  iterations 
for  convergence.  What  we  are  solving  for  are  corrections  to  our  a  priori  estimates  of  these  four 
quantities.  The  equation  of  condition  for  the  least-squares  solution,  which  is  computed  for  each 
observation,  is: 


fd He  df  dHc  dg  \  A  ,  fd He  df  dHc  dg  \  A  A 

\WWo  +  ~dXWo  )  A<^°  +  V  d<f>  d\0  +  ~dTd\^ )  AA° 

fdEcdf  HHc%\  fdEcdf  dEcdg\ 

v  d<t>  dC  +  dXdC'J  +  V  d<t>  dS  +  dX  dS  j 


where 


a  is  the  altitude  intercept  (a  =  Ho  —  He) , 

Ho  is  the  observed  altitude, 

He  is  the  computed  altitude, 

<f  is  the  estimated  latitude  for  the  time  of  observation, 

A  is  the  estimated  longitude  for  the  time  of  observation, 
<j)Q  is  the  estimated  latitude  for  the  time  of  the  fix, 

A0  is  the  estimated  longitude  for  the  time  of  the  fix, 

C  is  the  course  of  the  observer’s  vessel, 

S  is  the  speed  of  the  observer’s  vessel, 

/  is  the  rhumb-line  sailing  formula  for  latitude,  and 

g  is  the  rhumb-line  sailing  formula  for  longitude. 


(1) 


This  equation  expresses  the  altitude  intercept  a  in  terms  of  small  corrections  to  the  estimated 
latitude,  longitude,  course,  and  speed,  respectively  denoted  Af0,  AA0,  AC,  and  AS.  These  latter 
four  quantities  are  unknowns  to  be  determined,  and  once  their  values  are  known,  we  can  correct  our 
estimate  of  the  vessel’s  track.  The  corrections  to  latitude  and  longitude  have  0  subscripts  to  denote 
that  they  apply  to  the  specific  time,  t0,  for  which  the  fix  is  desired.  All  angles,  including  latitude 
and  longitude,  are  expressed  in  radians. 

The  quantities  within  parentheses  in  equation  (1)  are  the  coefficients  of  the  unknowns.  These 
coefficients,  along  with  the  altitude  intercept  a,  must  be  numerically  evaluated  for  each  observation. 
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To  do  that,  we  need  algebraic  expressions  for  the  partial  derivatives.  There  are  ten  partial  derivatives 
that  appear: 

dHc  dHc  df  df  df  df  dg  dg  dg  dg 

~W'  ~d\'  Wo'  W'  dC ’  dS’  Wo '  W'  dC ’  dS 

These  are  obtained  from  three  formulas:  the  equation  for  He  as  a  function  of  <f  and  A;  the  sailing 
formula  for  latitude,  represented  by  the  function  /;  and  the  sailing  formula  for  longitude,  represented 
by  the  function  g.  The  sailing  formulas  are  functions  that  together  provide  the  position  of  the  vessel 
as  a  function  of  time  f ;  they  involve  the  parameters  C ,  S,  fo,  and  A0  (the  sailing  formulas  are  derived 
in  [6]).  Algebraic  expressions  for  all  the  partials  needed  above  are  given  in  the  Appendix  to  [1], 

For  a  single-leg  solution,  the  mathematical  procedure  can  be  summarized  as  follows.  First,  the 
time  t0  for  the  fix  is  chosen.  Initial  estimates  of  the  values  of  f0,  A0,  C,  and  S  are  made  for  time 
t0.  The  sailing  formulas  /  and  g  then  allow  the  instantaneous  latitude,  <f>,  and  longitude,  A,  to  be 
computed  for  any  other  time  t.  For  each  observation,  therefore,  an  estimate  of  <f  and  A  is  made. 
The  sextant  altitude,  hs,  is  reduced  to  an  observed  altitude,  Ho,  in  the  usual  way.  An  almanac  or 
ephemeris  routine  is  then  used  to  provide  the  celestial  coordinates  (GHA  and  Dec)  of  the  observed 
body,  and  He  and  a  are  computed  for  the  observation.  Using  the  expressions  for  the  partials  given 
in  [1],  the  coefficients  of  the  unknowns  in  equation  (1)  can  be  computed.  To  include  the  observation 
in  the  least-squares  solution,  the  coefficients  are  added  as  a  new  row  in  the  design  matrix,  and  a 
is  added  to  the  column  vector  of  measurements.  After  all  observations  have  been  processed  in  this 
way,  the  least-squares  solution  is  computed,  yielding  values  for  Af0,  AA0,  AC,  and  AS.  These  are 
then  used  to  correct  the  original  estimates  of  f0,  A0,  C,  and  S  for  time  t0.  The  corrected  values 
of  (J)0  and  A0  define  the  fix.  The  whole  process  can  be  reiterated  if  necessary,  using  the  new  values 
of  (J)0,  Aq ,  C,  and  S  as  the  basis  for  the  computations.  Details,  including  a  numerical  example,  are 
given  in  [1],  For  a  description  of  least-squares  techniques,  see  any  text  on  scientific  data  analysis, 
for  example,  [7,8]. 

Equation  (1)  can  be  applied  beyond  the  case  of  celestial  navigation  if  we  simply  generalize  the 
meanings  of  several  variables.  Ho  and  He  can  be  regarded  as  the  observed  and  computed  values, 
respectively,  of  any  scalar  quantity  that  can  be  measured  while  underway  and  that  is  a  function  of 
the  observer’s  instantaneous  latitude  and  longitude.  For  example,  range  to  a  known  point  would 
be  such  a  quantity.  Then  a  =  Ho  —  He  becomes  simply  the  “observed  minus  computed”  residual  of 
this  quantity.  As  long  as  we  can  compute  a,  dUc/df,  and  dUc/dX  for  each  observation  (and  the 
geometry  of  the  observations  is  not  degenerate),  equation  (1)  applies.  The  development  presented 
below  does  not  depend  on  particular  meanings  of  Ho,  He,  a,  or  the  partial  derivatives  dUc/df  and 
dKc/dX. 


General  Considerations  for  Multiple-Leg  Tracks 

When  we  have  observations  spread  over  more  than  one  leg,  the  conditional  equation  gains  addi¬ 
tional  terms  corresponding  to  the  added  degrees  of  freedom  in  the  problem.  That  is,  the  geometry 
of  the  problem  becomes  more  complex  and  can  be  described  only  with  additional  parameters,  some 
of  which  we  must  solve  for.  We  will  consider  a  new  voyage  leg  to  begin  with  either  a  course  change 
or  a  speed  change  (or  both)  at  a  known  time.  We  model  these  changes  as  instantaneous,  and  we 
can  ignore  the  finite  turning  radius  of  the  vessel  if  (1)  no  observations  are  made  during  the  turn, 
and  (2)  the  turning  radius  is  short  compared  with  the  length  of  the  straight  portions  of  the  vessel’s 
track.  In  the  most  general  case,  each  leg  has  both  course  and  speed  independent  of  the  course  and 
speed  on  any  other  leg.  (This  allows  for  the  case  where  a  change  of  either  course  or  speed  alters 
the  interaction  of  wind  and  current  on  the  vessel  and  affects  the  other  parameter.)  Just  as  in  the 
single-leg  case,  we  assume  that,  initially,  we  know  the  vessel’s  track  reasonably  well  (errors  not 
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exceeding  a  few  degrees  in  position  or  course  and  a  few  knots  in  speed) ,  and  that  what  we  need  to 
solve  for  are  corrections  to  our  starting  estimates  of  the  relevant  parameters. 

To  cover  the  multiple-leg  case,  some  changes  to  our  mathematical  notation  are  necessary.  We 
will  use  subscripts  to  designate  specific  legs.  The  first  leg  will  be  denoted  by  subscript  1.  This  leg 
begins  at  time  tk,  when  the  vessel  is  at  position  ) ,  sailing  course  C\  at  speed  Sk.  The  second 

leg  begins  at  time  f2  when  the  vessel  is  at  position  (<(>2,A2);  the  new  course  and  speed  are  C 2  and 
SV  The  third  begins  at  time  f3,  and  so  on.  We  assume  that  we  know  the  times  tj  exactly,  and  that 
we  start  with  reasonable  estimates  for  the  other  quantities. 

We  have  used  the  symbols  /  and  g  to  represent  the  sailing  formulas  for  rhumb-line  tracks,  which 
respectively  provide  the  latitude  and  longitude  of  the  vessel  as  a  function  of  time.  For  the  multiple- 
leg  case,  we  must  use  the  symbols  and  gk  to  represent  the  sailing  formulas  to  be  used  for  leg  1,  /2 
and  g2  to  represent  the  sailing  formulas  to  be  used  for  leg  2,  etc.  The  algebraic  form  of  the  sailing 
formulas  is,  of  course,  the  same  for  all  legs;  it  is  just  the  parameter  values  used  in  the  formulas 
that  are  different  for  each  leg.  We  will  continue  to  use  the  symbols  /  and  g  (without  subscripts)  to 
represent  the  generic  form  of  the  sailing  formulas.  The  values  of  five  parameters  in  these  formulas 
effectively  define  the  functions  fj  and  gj  for  each  leg  j.  The  parameters  are  the  time  tj  at  which  the 
leg  begins,  the  latitude  epj  and  longitude  A  j  of  the  vessel  at  time  tj,  the  course  Cj,  and  the  speed 
Sj.  For  each  leg  j  we  have 


H*)  =fj  =  f(t~  ^  XjiCjiSj) 
A  (t)  gj  g  (t  tj ,  <pj ,  A  j ,  Cj ,  Sj ) 


(2) 


The  initial  point  of  the  leg  (position  ((pj,Xj)  at  time  tj)  thus  serves  as  a  reference  point  for  the  sailing 
formulas  for  the  leg.  We  have  made  the  reference  point  of  each  leg  to  be  the  initial  point,  but  that 
is  a  requirement  only  for  the  second  and  succeeding  legs  in  a  multiple-leg  case;  for  the  first  (or  only) 
leg,  the  reference  point  can  be  anywhere  along  the  leg. 

It  should  be  noted  that  the  functional  dependence  of  the  two  sailing  formulas  on  the  parameters 
shown  in  equation  (2)  is  implicit;  the  set  of  parameters  used  explicitly  in  the  two  formulas  is 
somewhat  different.  As  presented  in  [1]  and  [6],  both  formulas  contain  auxiliary  variables  that  need 
to  be  expanded  before  the  dependence  of  these  formulas  on  the  parameter  set  shown  above  becomes 
evident. 

The  naive  approach  to  the  multiple-leg  problem  would  be  to  simply  take  equation  (1)  and 
generalize  it  to  cover  N  legs,  using  the  notation  just  introduced: 


N 

3  =  1 


+ 


cddc  dfk  cddc  dgk ' 
d(j>  d(j>j  <9  A  dcpj, 

cddc  dfk  i  cddc  dgk 
d<j>  dC3  +  d\  dCj 


A  <t>j  + 

ACj  + 


<9Hc  dfk  <9Hc  dgk\  A  A 
d<t>  d\j  dX  dXj)  3 

3dEcdfk  <9Hc  dgk\  A 
v  d<t>  dS3  +  dX  dSj)  3 


(3) 


where  the  subscript  k  denotes  the  leg  on  which  the  observation  was  taken,  a^k  represents  the  altitude 
intercept  for  the  zth  observation  on  leg  k ,  and  the  summation  is  over  all  N  legs  in  the  problem.  There 
are  4N  free  parameters  to  be  determined:  A  cpj,  AA  j,  ACj,  and  A  Sj  for  j  =  1  through  N.  The  partial 
derivatives  of  the  sailing  formulas  that  appear  in  equation  (3)  contain  subscripts.  The  subscript  in 
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the  “numerator”,  k,  represents  the  leg  number  of  the  observation  being  processed  and  the  subscript 
in  the  “denominator”,  j,  represents  the  leg  number  of  a  particular  unknown. 

In  this  construction  of  the  problem,  each  leg  is  independent  of  the  others,  and  there  is  no 
requirement  that  the  legs  be  continuous.  With  independent  legs,  only  observations  taken  on  a 
particular  leg  contribute  to  the  parameters  for  that  leg.  That  is  because  the  sailing  formulas  for  a 
particular  leg  depend  only  on  parameters  relating  to  that  leg,  and  therefore  the  partial  derivatives 
of  the  sailing  formulas  are  zero  if  j  /  k.  Therefore,  in  equation  (3),  as  j  cycles  from  1  to  N,  the 
coefficients  of  all  the  unknowns  are  zero  except  when  j  =  k,  so  that  an  observation  has  no  effect 
except  on  parameters  relating  to  its  own  leg.  When  j  =  k,  the  partials  are  evaluated  from  the  single¬ 
leg  formulas  given  in  [1],  What  we  have  done,  therefore,  is  set  up  a  single  least-squares  solution 
for  a  problem  that  could  be  handled  as  N  individual  least-squares  solutions,  one  for  each  leg.  The 
values  of  the  unknowns  that  we  solve  for  will  be  the  same  regardless  of  which  approach  is  taken. 

Of  course,  the  idea  that  a  multiple-leg  voyage  consists  of  a  series  of  independent  legs  does  not 
correspond  to  the  real  problem.  The  courses  and  speeds  may  be  independent,  but  each  leg  must  be 
continuous  with  its  neighbors.  Because  we  are  modeling  each  leg  change  as  instantaneous,  continuity 
is  required  only  in  position,  not  in  motion.  The  requirement  that  adjacent  legs  be  continuous  defines 
a  set  of  constraints  on  the  solution,  and  constraints  always  reduce  the  number  of  free  parameters  in 
a  problem.  Formal  methods  have  been  developed  to  deal  with  constrained  least-squares  problems, 
but  the  multiple-leg  problem  can  be  attacked  more  easily.  In  the  next  section,  we  will  abandon 
equation  (3)  and  take  another  approach  to  the  problem. 


Algorithm  for  Multiple-Leg  Tracks 

A  viable  solution  to  the  real  multiple-leg  problem  must  explicitly  or  implicitly  deal  with  the 
requirement  of  continuity  at  the  leg  intersections.  The  mathematical  condition  is  that,  for  all  legs 
beyond  the  first,  the  point  at  the  start  of  leg  j  must  be  identical  to  that  at  the  end  of  leg  j  —  1.  At 
the  start  of  leg  j  >  1,  the  vessel’s  coordinates  are  therefore  not  free  parameters.  Only  on 

the  first  leg  are  the  coordinates  of  the  reference  point  free  parameters.  Given  the  location  of  that 
reference  point,  the  time  of  each  course  or  speed  change,  and  the  course  and  speed  values  on  each 
leg,  the  vessel’s  track  is  completely  determined.  Since  all  the  times  are  known,  for  N  legs,  there  are 
2N+2  free  parameters  in  the  problem  (rather  than  4N,  as  in  equation  (3)  above).  We  start  with 
estimates  of  their  values  and  seek  to  correct  those  estimates.  The  two  quantities  (pi  and  Ai,  which 
define  the  position  of  the  vessel  at  the  reference  point  of  the  first  leg,  represent  the  fix.  The  other 
free  parameters  of  the  problem  are  the  course  and  speed  along  each  leg.  The  conditional  equation 
for  this  case  must  therefore  take  the  form 


N 

+£ 

j=i 


<9Hc  dfk  <9Hc  dgk\  A  x 

df  dfa  d\  dfj  n 

<9Hc  dfk  <9Hc  dgk\ 

d<t>  dCj  d\  dCj )  3 


f  <9Hc  dfk  <9Hc  dgk  \  .  , 

/  <9Hc  dfk  <9Hc  dgk  \  A  c, 
V  d<t>  dS3  +  d\  dSj)  3 


(4) 


where  the  notation  is  the  same  as  for  equation  (3) .  This  equation  must  be  evaluated  numerically  for 
each  observation,  and  the  observations  are  taken  in  chronological  order  so  that  all  the  observations 
from  a  given  voyage  leg  k  are  processed  before  any  observations  from  the  next  leg. 

In  equation  (4),  as  in  equations  (1)  and  (3),  the  coefficients  of  the  unknowns  are  composed  of 
linear  combinations  of  partial  derivatives.  The  issue  is  how  to  evaluate  these  partial  derivatives  as 
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the  observations  from  each  leg  are  processed.  These  partial  derivatives  are  of  two  types.  The  partials 
of  the  observable  He  with  respect  to  the  instantaneous  position  of  the  vessel  (dUc/dcp  and  dUc/dX) 
are  straightforward  to  evaluate;  if  He  is  celestial  altitude,  the  expressions  can  be  taken  from  [1],  The 
partials  of  the  sailing  formulas  /*.  and  gj~  are  more  difficult  to  deal  with,  because  they  describe  the 
relationship  between  the  ship’s  instantaneous  position  and  all  the  free  parameters  in  the  problem 
((pi  and  Xi  as  well  as  Cj  and  Sj  for  all  j). 

The  sailing  formula  partials  contain  subscripts:  k  in  the  “numerator”  represents  the  leg  number 
of  the  current  observation  and  j  in  the  “denominator”  represents  the  leg  number  of  a  particular 
unknown.  (In  the  top  line  of  equation  (4),  j  =  1.)  Within  each  evaluation  of  equation  (4)  (each 
observation),  k  remains  fixed  but  j  runs  over  all  the  legs  in  the  problem.  In  evaluating  the  sailing- 
formula  partials,  we  therefore  have  three  cases  to  consider:  j  <  k,  j  =  k,  and  j  >  k.  Keeping  in 
mind  the  meanings  of  j  and  k,  we  recognize  that  for  j  >  k,  the  sailing  formula  partials  must  be 
zero,  reflecting  the  fact  that  the  ship’s  current  computed  position  cannot  depend  on  the  parameters 
of  future  legs.  For  j  =  k,  the  partials  must  be  the  same  as  in  the  single-leg  case,  and  the  expressions 
from  [1]  can  be  directly  used.  The  j  <  k  case  is  the  most  complicated,  because  these  partials 
must  reflect  the  fact  that  the  ship’s  computed  position  does  depend  on  the  parameters  of  past  legs, 
propagated  through  the  constraints  at  the  leg  intersections.  Most  of  Appendix  A  deals  with  this 
case. 

Thus,  the  net  effect  of  the  partials  of  the  sailing  formulas  on  equation  (4)  is  that  each  obser¬ 
vation  contributes  to  determining  not  just  the  course  and  speed  on  its  own  leg  k  but  also  those 
of  all  preceding  legs.  However,  the  future  is  indeterminate,  and  observations  cannot  provide  any 
information  on  the  unknowns  relating  to  legs  that  have  not  yet  occurred  (they  may  not  occur) .  All 
observations  contribute  to  determining  the  coordinates  of  the  reference  point  on  the  first  leg,  which 
serves  as  the  fix  for  the  multiple-leg  problem. 

In  Appendix  A,  an  algorithm  for  the  evaluation  of  the  sailing  formulas  partials  that  appear  in 
equation  (4)  is  presented.  The  algorithm  is  quite  straightforward  in  that  all  that  is  done  for  the  j  <  k 
case  is  to  implicitly  fold  the  continuity  constraints  into  the  sailing  formulas,  the  functions  /*.  and  gj~ 
for  leg  k.  Then,  new  expressions  for  the  partial  derivatives  of  /*.  and  gj~  are  immediately  obtained 
using  the  chain  rule.  These  expressions  amount  to  recurrence  relations.  Partials  computed  for  each 
leg  are  built  on  partials  computed  for  previous  legs  (a  process  that  follows  logically  from  the  nature 
of  the  problem).  Some  care  with  mathematical  bookkeeping  is  necessary.  However,  the  advantage 
of  this  scheme  is  that,  if  the  observations  are  processed  in  chronological  order,  no  calculations  are 
required  beyond  those  needed  for  the  single-leg  case,  except  trivial  multiplications  and  additions. 
No  changes  are  required  in  the  least-squares  procedure  itself.  Appendix  A  describes  the  procedure 
in  detail.  Because  it  represents  a  modification  to  the  single-leg  case,  the  algorithm  is  not  difficult  to 
implement  in  computer  code. 

In  the  development  of  the  multiple-leg  case  presented  here,  the  fix  is  always  on  the  first  leg.  In 
practice,  what  is  of  interest  is  not  really  the  vessel’s  position  on  the  first  leg,  but  the  last,  which  is 
the  leg  that  the  vessel  is  currently  sailing.  It  may  seem  as  if  the  problem  should  be  reversed.  In 
fact,  we  can  do  just  that — the  problem  is  symmetric  in  time  and  the  vessel  can,  mathematically, 
be  sailed  backwards  as  well  as  forwards.  Therefore,  what  we  have  been  calling  the  “first  leg”  can 
actually  be  the  last  and  what  we  called  the  “initial  point”  of  each  leg  can  actually  be  the  final 
point.  The  problem  has  been  described  in  a  manner  that  keeps  the  development  and  terminology 
as  straightforward  as  possible.  This  allows  us  to  visualize  the  vessel  as  moving  forward  in  time,  to 
discuss  causes  and  effects  in  a  natural  order,  and  to  refer  to  leg  j  +  1  as  the  “subsequent”  or  “next” 
leg.  However,  when  the  algorithm  presented  here  is  used  for  real  problems,  the  procedure  can  be 
implemented  in  a  time-reversed  fashion  to  provide  information  of  practical  use. 
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Sample  Calculation 


A  numerical  example  using  this  algorithm  is  presented  below.  A  set  of  12  artificial  sextant 
observations  (hs  values),  spanning  3  legs  of  a  hypothetical  voyage  (4  observations  per  leg),  was 
reduced  using  the  algorithm  described  in  this  paper.  The  solution  yielded  a  position  fix  for  a  single 
specified  time  in  the  final  leg  and  course  and  speed  estimates  for  all  legs.  The  artificial  observations 
were  generated  by  a  computer  program  that  numerically  integrated  a  hypothetical  ship’s  course, 
consisting  of  3  rhumb-line  tracks,  as  shown  in  Table  1  (north  latitudes  and  east  longitudes  are 
positive) . 


Table  1 — Ship’s  Track  Used  in  Sample  Calculation 


Start  of  Leg  A 

Start  of  Leg  B 

Start  of  Leg  C 

Date: 

30  Aug  1996 

Date: 

31  Aug  1996 

Date: 

31  Aug  1996 

Time: 

14:00:00  UTC 

Time: 

00:18:30  UTC 

Time: 

14:06:26  UTC 

Latitude: 

+40.000° 

Latitude: 

+43.275° 

Latitude: 

+46.446° 

Longitude: 

-50.000° 

Longitude: 

-47.479° 

Longitude: 

-43.736° 

Course: 

30.00° 

Course: 

40.00° 

Course: 

70.00° 

Speed: 

22.00  kn 

Speed: 

18.00  kn 

Speed: 

15.00  kn 

The  data  in  Table  1,  and  any  positions  derived  from  them,  constitute  “truth”  for  this  case.  The 
latitude  and  longitude  listed  for  the  beginning  of  legs  B  and  C  were  computed  from  the  data  for 
the  previous  leg;  that  is,  the  legs  meet  at  these  points  and  the  changes  to  course  and  speed  were 
assumed  to  be  instantaneous.  Astronomical  positions  of  the  celestial  bodies  chosen  were  computed 
that  are  consistent  with  those  in  the  1996  Nautical  Almanac  and  transformed  to  local  altitudes. 
Atmospheric  refraction  was  included  but  the  height  of  eye  was  considered  to  be  zero,  and  the  time 
scale  UT1  was  assumed  to  be  identical  to  UTC  (that  is,  DUT1=UT1-UTC=0).  Random  errors  with 
a  0.7  arcmin  standard  deviation  were  added  to  the  calculated  hs  values.  The  artificial  observations, 
for  30-31  August  1996  (UTC),  are  given  in  Table  2. 


Table  2 — Artificial  Observations  for  Sample  Calculation 


Leg  A  —  30  Aug 

Leg  B 

—  31  Aug 

Leg  C  —  31  Aug 

UTC 

Object 

hs 

UTC 

Object 

hs 

UTC 

Object 

hs 

14:58:54.3 

Sun 

57.872 

00:36:29.7 

Moon 

17.659 

22:01:35.5 

Alpheratz 

19.411 

18:25:04.2 

Sun 

36.890 

05:42:28.2 

Moon 

48.486 

22:23:11.4 

Rasalhague 

54.089 

22:50:08.4 

Arcturus 

34.624 

06:54:27.9 

Markab 

37.710 

22:44:47.3 

Kochab 

53.487 

23:04:51.9 

Jupiter 

23.795 

11:42:26.5 

Sun 

33.730 

23:06:23.2 

Alpheratz 

30.225 

All  of  the  Sun  and  Moon  observations  apply  to  the  lower  limb.  The  star  and  planet  observations 
represent  twilight  observations;  the  two  Moon  observations  represent  night  observations.  The  Moon’s 
phase  was  0.92  (age  17  days). 

Once  computed  and  stored  in  a  hie,  these  artificial  observations  were  used  as  input  to  a  separate 
computer  program  that  implemented  the  algorithm  described  in  this  paper.  This  second  program 
was  given  the  incorrect  data  on  the  ship’s  position  and  motion  listed  in  Table  3. 


7 


Table  3 — Incorrect  Data  on  Ship’s  Track 
Supplied  to  Multiple-Leg  Algorithm 


Start  of  Leg  A 

Start  of  Leg  B 

Start  of  Leg  C 

Date: 

30  Aug  1996 

Date: 

31  Aug  1996 

Date: 

31  Aug  1996 

Time: 

14:00:00  UTC 

Time: 

00:18:30  UTC 

Time: 

14:06:26  UTC 

Latitude: 

+41.700° 

Latitude: 

+44.469° 

Latitude: 

+47.159° 

Longitude: 

-48.000° 

Longitude: 

-45.639° 

Longitude: 

-42.050° 

Course: 

32.00° 

Course: 

43.00° 

Course: 

67.50° 

Speed: 

19.00  kn 

Speed: 

16.00  kn 

Speed: 

12.00  kn 

The  dates  and  times  shown  in  Table  3  are  the  same  as  in  Table  1,  but  the  course  and  speed  values, 
as  well  as  the  position  of  the  first  point  of  leg  A,  are  considerably  different.  Just  as  in  Table  1,  the 
latitude  and  longitude  listed  in  Table  3  for  the  beginning  of  legs  B  and  C  are  computed  from  the 
data  listed  for  the  previous  leg.  Clearly,  this  track  (Table  3)  is  a  rather  poor  approximation  to  the 
true  track  (Table  1).  The  test  of  the  algorithm  is  whether,  given  only  the  12  observations  (Table  2), 
it  can  determine  the  true  track  of  the  vessel,  limited  only  by  the  error  in  the  observations.  The  date 
and  time  for  the  desired  fix  were  chosen  to  be  1  Sep  1996  at  07:30:00  UTC,  which  is  on  leg  C,  8.4 
hours  after  the  last  observation.  In  addition  to  the  latitude  and  longitude  for  that  time,  course  and 
speed  for  all  3  legs  were  to  be  determined.  Based  on  the  incorrect  data  on  the  ship’s  track  it  was 
given,  the  second  program  computed  the  estimated  latitude  and  longitude  for  the  time  of  the  fix  to 
be  (+48.490°,  —37.280°).  It  also  computed  the  data  in  Table  4  for  the  12  observations. 


Table  4 — Estimates  of  Sight-Reduction  Quantities  for  Observations 


Leg 

UTC 

Est  (j) 

Est  A 

Object 

Ho 

He 

a 

0 

0 

0 

0 

/ 

A 

14:58:54.3 

+41.964 

-47.780 

Sun 

57.862 

56.407 

+87.3 

A 

18:25:04.2 

+42.887 

-47.001 

Sun 

36.868 

34.944 

+115.5 

A 

22:50:08.4 

+44.073 

-45.983 

Arcturus 

34.600 

33.126 

+88.4 

A 

23:04:51.9 

+44.139 

-45.926 

Jupiter 

23.758 

22.465 

+77.6 

B 

00:36:29.7 

+44.527 

-45.563 

Moon 

17.607 

18.608 

-60.0 

B 

05:42:28.2 

+45.522 

-44.255 

Moon 

48.472 

47.192 

+76.8 

B 

06:54:27.9 

+45.756 

-43.944 

Markab 

37.688 

36.198 

+89.4 

B 

11:42:26.5 

+46.692 

-42.687 

Sun 

33.706 

34.379 

-40.4 

C 

22:01:35.5 

+47.765 

-39.893 

Alpheratz 

19.364 

20.332 

-58.1 

C 

22:23:11.4 

+47.793 

-39.795 

Rasalhague 

54.077 

53.217 

+51.5 

c 

22:44:47.3 

+47.820 

-39.696 

Kochab 

53.474 

53.770 

-17.8 

c 

23:06:23.2 

+47.848 

-39.597 

Alpheratz 

30.197 

31.023 

-49.6 

In  Table  4,  the  altitude  intercept,  a,  is  Ho  —  He  and  positive  values  indicate  “toward.”  The  Ho  and 
He  values  are  topocentric  and,  for  the  Sun  and  Moon,  they  refer  to  the  lower  limb.  The  large  a 
values  are  a  result  of  the  poor  data  on  the  ship’s  track  that  the  program  was  given,  that  we  hope  to 
correct. 

As  described  previously,  the  algorithm  described  in  this  paper  will  usually  be  implemented  in 
a  time-reversed  sense.  The  program  followed  this  strategy,  working  backwards  through  the  obser¬ 
vations  and  effectively  treating  leg  C,  containing  the  fix,  as  the  “first”  leg.  In  recognition  of  this 
strategy,  and  in  order  to  maintain  consistency  with  equation  (4) ,  the  following  notation  will  therefore 
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be  used: 


Leg  A  — >  Leg  3:  Course  C3,  Speed  S3 
Leg  B  — >  Leg  2:  Course  C2,  Speed  S2 

Leg  C  — >  Leg  1:  Course  C\,  Speed  S 1,  Fix  Coordinates  (pi,  Xi 

As  it  processed  each  observation,  the  program  computed  the  values  of  the  coefficients  of  all  the 
unknowns  in  equation  (4).  It  sent  these  coefficient  values,  along  with  the  value  of  a,  to  a  standard 
least-squares  routine.  All  the  observations  were  given  equal  weight.  After  all  12  observations  had 
been  processed  in  this  way,  the  least-squares  solution  was  formed.  The  results  are  given  in  Table  5, 
which  lists  the  corrections  to  the  parameters  describing  the  ship’s  track. 


Table  5 — Least-Squares  Solution 


Parameter 

No.  Obs. 

Value 

A</>! 

12 

-0.579  ±0.041 

O 

AAj 

12 

-0.517  ±0.045 

0 

ACi 

12 

±2.76 

±0.85 

0 

ASi 

12 

±2.77 

±0.20 

kn 

ac2 

8 

-3.01 

±0.38 

O 

as2 

8 

±2.25 

±0.21 

kn 

AC's 

4 

-2.30 

±0.48 

O 

as3 

4 

±2.87 

±0.20 

kn 

The  second  column  in  Table  5  indicates  the  number  of  observations  contributing  to  each  parameter. 
Table  5  shows  rather  large  corrections  to  the  parameter  values,  as  we  expect  from  the  crudeness  of 
the  track  data  that  the  program  used.  The  mean  error  of  unit  weight  of  the  solution  was  computed 
to  be  0.8  arcmin.  If  we  add  the  parameter  values  from  Table  5  to  the  estimated  (incorrect)  data  on 
the  ship’s  track,  we  obtain  the  final  results  of  the  entire  process,  summarized  in  Table  6. 

Table  6 — Results  of  Solution  for  Ship’s  Track 


Truth 

Solution 

Error 

i>\ 

±47.932 

±47.911 

-0.021 

Ai 

-37.745 

-37.797 

-0.052 

Cr 

70.00 

70.26 

±0.26 

Si 

15.00 

14.77 

-0.23 

C2 

40.00 

39.99 

-0.01 

s2 

18.00 

18.25 

±0.25 

C  3 

30.00 

29.70 

-0.30 

53 

22.00 

21.87 

-0.13 

The  units  used  in  Table  6  are  degrees  and  knots,  as  in  the  previous  tables.  “Error”  simply  refers  to 
the  difference  solution-truth.  Obviously  a  great  improvement  in  the  data  on  the  ship’s  track  was 
made.  Of  course,  in  a  real  situation,  “Truth”  would  not  be  known — the  question  at  this  point  would 
be  whether  another  iteration  of  entire  procedure  would  be  advisable,  using  the  updated  values  of 
the  parameters  describing  the  ship’s  track  given  in  Table  6.  In  fact,  the  program  did  perform  two 
iterations  of  the  algorithm,  and  the  results  given  in  Tables  5  and  6  reflect  the  sum  of  the  corrections 
obtained  from  the  two  iterations.  If  the  solution  is  reliable,  the  iterations  will  converge  rapidly;  that 
is,  the  parameter  corrections  from  a  second  or  third  iteration  will  all  be  much  smaller  than  those 
from  the  previous  iteration.  That  was  the  case  here. 
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The  algorithm  did  recover  the  ship’s  true  track,  although  the  fix  coordinates  are  perhaps  not  as 
accurate  as  we  might  desire.  The  error  in  the  fix  position  is  2.5  nmi,  even  though  the  observations 
had  a  scatter  of  0.7  arcmin.  This  example  is  adversely  affected  by  the  small  number  of  observations 
involved — there  are  only  times  as  many  observations  as  there  are  parameters,  leaving  only  4 
degrees  of  freedom.  The  parameter  correlation  matrix  from  the  solution  showed  a  0.96  correlation 
between  and  A C^,  indicating  that  the  geometry  of  the  observations  does  not  allow  a  clean 
separation  of  the  effects  of  these  errors.  Such  problems  are  not  uncommon  in  cases  with  so  few 
observations.  Generally,  at  least  a  half-dozen  observations  on  each  leg,  well  spread  out  in  time, 
should  be  used.  In  one  scenario  tested,  adding  just  4  more  observations  to  this  example  reduced  the 
error  in  the  fix  to  less  than  1  nmi.  The  extra  observations  not  only  doubled  the  number  of  degrees 
of  freedom  in  the  problem  but  also  improved  its  geometry. 


In  any  event,  the  multiple-leg  algorithm  described  in  this  paper  may  provide  the  best  answers 
attainable  for  the  12-observation  case  we  have  been  considering.  For  example,  if  we  try  something 
more  conventional,  and  process  only  the  observations  in  leg  1  (leg  C),  we  obtain  the  results  in 
Table  7. 


Table  7 — Results  of  Solution  for  Leg  1  Only 


Truth  Solution 

+47.932  +47.847 

Aj  -37.745  -38.411 


Error 

-0.085 

-0.666 


fa  +47.932 

Aj  -37.745 

C\  70.00 

Sj  15.00 


+47.720  -0.212 
-37.500  +0.245 
76.45  +6.45 

15.76  +0.76 


Two  solutions  for  leg  1  (leg  C)  are  given  in  Table  7:  above  the  line,  only  the  latitude  and  longitude 
of  the  fix  have  been  solved  for;  below  the  line,  course  and  speed  values  have  also  been  determined 
(the  second  solution  has  no  degrees  of  freedom) .  Clearly,  neither  is  satisfactory — the  errors  in  the  fix 
position  are  27.6  and  16.1  nmi,  respectively.  The  difficulty  in  treating  leg  1  alone  is  that  the  obser¬ 
vations  represent  a  group  of  sights  taken  within  a  1-hour  period  during  evening  twilight.  Since  the 
course  and  speed  are  not  well  known,  and  cannot  be  accurately  determined  using  these  observations, 
the  extrapolation  to  the  time  of  the  fix  fails.  In  this  case,  the  use  of  observations  from  the  previous 
legs  can  add  substantial  information  to  the  process,  providing  a  reasonable  determination  of  the 
course  and  speed  values  for  leg  1  (the  final  leg) ,  and  allowing  a  relatively  accurate  extrapolation  of 
the  ship’s  position  forward. 


What  would  happen  to  our  3-leg,  12-observation  case  if  the  observations  were  substantially 
better  (or  worse)  than  those  listed  in  Table  2?  We  can  easily  change  the  random  scatter  in  the 
artificial  observations  and  run  through  the  entire  procedure  again.  The  results  of  this  experiment 
are  summarized  in  Table  8,  the  body  of  which  shows  the  error  of  each  parameter  (solution-truth), 
in  the  units  we  have  been  using,  as  a  function  of  observational  scatter. 
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Table  8 — Parameter  Accuracy  vs.  Observational  Accuracy 


—  Observational  Uncertainty  — 


0 

±0.05' 

±0.2' 

±0.7' 

±1.5' 

i>\ 

+0.001 

-0.001 

-0.006 

-0.021 

-0.047 

Ai 

+0.001 

-0.003 

-0.015 

-0.052 

-0.113 

Cr 

-0.02 

0.00 

+0.06 

+0.26 

+0.60 

Si 

0.00 

-0.02 

-0.06 

-0.23 

-0.49 

c2 

0.00 

0.00 

0.00 

-0.01 

-0.02 

s2 

0.00 

+0.01 

+0.07 

+0.25 

+0.54 

c3 

0.00 

-0.02 

-0.08 

-0.30 

-0.64 

s3 

0.00 

-0.01 

-0.04 

-0.13 

-0.27 

There  are  no  surprises  in  Table  8;  it  simply  indicates  that  the  errors  in  the  parameters  solved  for 
are  directly  proportional  to  the  scatter  in  the  observations. 


Use  of  the  Multiple-Leg  Algorithm 

The  algorithm  outlined  above  and  detailed  in  Appendix  A  has  been  implemented  in  software 
for  celestial  navigation.  It  was  produced  by  the  U.S.  Naval  Observatory  for  Navy  shipboard  use  [9]. 
This  software  is  currently  operational  in  the  fleet.  The  algorithm  has  been  extensively  tested  using 
artificial  observations  generated  by  independent  programs,  where  “truth”  can  be  known  absolutely 
and  controlled.  The  algorithm  is  accurate  and  remarkably  robust  across  a  wide  variety  of  voyage 
track  scenarios. 

The  algorithm  cannot,  however,  be  applied  blindly,  and  the  distribution  of  the  observations  in 
time  and  azimuth  should  be  considered  before  proceeding.  As  the  sample  calculation  above  shows, 
it  is  advisable  for  there  to  be  at  least  a  half-dozen  observations  on  each  leg.  However,  valid  solutions 
can  sometimes  be  formed  with  fewer  observations,  or  even  when  one  or  more  legs  are  “empty.”  The 
example  above  also  shows  the  utility  of  the  algorithm  in  establishing  a  mathematical  linkage  between 
the  legs  that  allows  observations  from  subsequent  (previous)  legs  to  contribute  to  the  solution  on  the 
first  (last)  leg — which  includes  the  adjustment  to  the  fix  coordinates.  It  is  difficult  to  provide  hard 
and  fast  a  priori  rules  for  the  algorithm’s  application  that  apply  in  all  cases.  As  in  any  least-squares 
procedure,  the  parameter  correlation  matrix  (obtained  from  the  variance-covariance  matrix)  from 
the  solution  should  be  examined  to  determine  whether  the  solution  is  reliable.  When  that  is  not 
the  case,  it  may  be  inadvisable  to  use  the  algorithm  in  its  fullest  generality  as  presented  above; 
course  and  speed  corrections  for  certain  legs  may  be  excluded.  In  such  cases,  in  the  summation 
term  of  equation  (4),  the  j  for  any  leg  can  simply  be  skipped  over  and  the  corresponding  A Cj 
and  ASj  not  included  in  the  solution.  If  the  distribution  of  observations  is  very  sparse,  the  entire 
summation  term  of  equation  (4)  can  be  omitted,  and  the  solution  confined  to  A<j)\  and  AAi,  that  is, 
only  corrections  to  the  fix  coordinates  are  computed.  Clearly,  however,  the  algorithm  is  most  useful, 
and  reliable,  when  there  is  a  good  distribution  of  observations  on  each  leg.  In  the  case  of  celestial 
navigation,  an  automated  observing  system  might  be  required  to  provide  the  number  and  accuracy 
of  the  observations  needed. 

The  algorithm  is  quite  general.  As  previously  noted,  equation  (4)  is  not  limited  to  celestial 
observations — He  can  represent  any  observable  which  is  a  known  function  of  latitude  and  longitude. 
The  only  requirement  is  that  the  functions  /  and  g  provide  the  coordinates  of  the  vessel  as  a  function 
of  time  and  the  latitude,  longitude,  course,  and  speed  at  the  beginning  of  the  leg. 
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The  algorithm  might  be  useful  even  when  a  ship  is  nominally  holding  a  fixed  course  and  speed. 
If  significant  changes  in  current  or  wind  occur,  a  single  rhumb  line  may  not  adequately  model  the 
ship’s  track  over  bottom.  In  such  cases,  the  ship’s  motion  might  be  approximated  by  a  series  of 
rhumb-line  legs,  and  the  algorithm  given  in  this  paper  would  then  be  applicable. 

Finally,  in  [1],  the  navigation  problem  for  vessels  on  the  surface  of  the  Earth  was  described 
as  being  analogous  to  the  orbit  determination  problem  for  bodies  in  the  solar  system.  That  paper 
derived  the  conditional  equation  for  the  single-leg  navigation  problem,  reproduced  here  as  equa¬ 
tion  (1),  from  considerations  similar  to  those  used  by  astronomers  for  the  orbit  problem.  It  may 
seem  as  if  the  multiple-leg  navigation  problem  addressed  in  this  paper  has  no  analog  in  orbit  theory, 
but  that  is  not  so.  Spacecraft  trajectories  are  adjusted  by  rocket  brings,  and  the  orbits  of  asteroids, 
comets,  and  spacecraft  are  perturbed  by  close  approaches  to  massive  bodies.  In  most  of  these  cases 
the  agent  of  change  (a  rocket  burn  or  close  gravitational  encounter)  is  a  transient  phenomenon.  An 
approximate  treatment  of  such  a  problem,  adequate  for  many  purposes,  would  be  to  consider  the 
“before”  and  “after”  trajectories  to  meet  at  a  point  in  space  where  the  motion  of  the  body  instanta¬ 
neously  changes.  Thus,  even  when  orbits  in  the  solar  system  are  considered,  we  have  situations  that 
are  quite  similar  to  the  multiple-leg  navigation  case.  The  strategy  for  dealing  with  the  multiple-leg 
navigation  problem  described  in  this  paper  could  be  generalized  to  the  orbit-determination  problem 
as  well. 


Conclusion 

The  running  fix — the  determination  of  position  at  sea  using  observations  taken  over  extended 
time  spans — is  a  common  navigational  problem.  Traditional  chart-based  techniques  advance  or 
retire  lines  of  position  to  a  common  time.  The  technique  works  best  when  the  vessel’s  motion  is 
constant,  but  it  is  also  applied  across  changes  to  course  or  speed,  that  is,  to  multiple-leg  tracks.  A 
commonly  used  mathematical  development  of  celestial  navigation  mimics  the  chart-based  scheme. 
The  weakness  of  this  approach  to  the  running-fix  problem  is  that  it  is  based  on  the  assumption  that 
the  difference  between  a  vessel’s  estimated  and  true  positions  remains  constant  in  two  dimensions 
as  the  ship  moves.  The  notion  that  a  ship’s  positional  error  is  constant  is  unlikely  to  be  true,  even 
approximately,  except  over  very  short  periods,  and  any  change  to  course  or  speed  only  complicates 
the  situation. 

An  algorithm  has  been  presented  in  this  paper  that  allows  celestial  or  other  observations  made 
over  several  voyage  legs  to  be  combined  into  a  single  solution  for  all  relevant  sailing  parameters:  the 
latitude  and  longitude  of  one  point  on  the  track  and  the  course  and  speed  along  each  leg.  In  this 
algorithm,  the  error  in  position  is  assumed  to  change,  and  its  rate  of  change  is  assumed  to  be  different 
along  each  leg.  The  algorithm  is  conceptually  simple,  uses  ordinary  least-squares  procedures,  makes 
efficient  use  of  the  available  observations,  and  is  relatively  easy  to  implement  in  computer  code. 
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Appendix  A 

Evaluating  the  Partial  Derivatives  of  the  Sailing  Formulas  in  the  Multiple-Leg  Case 


Equation  (4)  in  the  main  part  of  the  paper  is  the  conditional  equation  for  a  least-squares 
solution  for  the  unknowns  of  the  multiple-leg  navigation  problem.  To  use  it,  we  need  to  be  able  to 
numericahy  evaluate  the  quantities  within  the  parentheses  in  the  equation,  which  are  the  coefficients 
of  the  unknowns.  Equation  (4)  is  evaluated  for  each  observation.  It  is  advantageous  to  take  the 
observations  in  chronological  order,  so  that  all  observations  on  a  given  leg  are  processed  before  any 
observations  on  the  next  leg. 

The  coefficients  of  the  unknowns  in  equation  (4)  consist  of  linear  combinations  of  various  partial 
derivatives.  As  mentioned  in  the  main  text,  the  partials  that  appear  are  either  partials  of  the 
observable  He  or  partials  of  the  sailing  formulas  /  and  g.  For  celestial  observations,  expressions 
for  the  partial  derivatives  dUc/dcp  and  dUc/dX  are  given  in  [1]  and  will  not  be  discussed  further 
here.  For  the  multiple-leg  case,  the  quantities  that  need  special  attention  are  the  partial  derivatives 
of  the  sailing  formulas,  for  example,  dfk/dfii  or  dgk/dCj.  These  appear  in  equation  (4)  in  the 
generic  form  dyijdxj,  where  y  represents  either  the  function  /  or  g,  and  x  represents  any  of  the  four 
parameters  <f>,  A,  C,  or  S.  We  will  refer  to  these  partial  derivatives  simply  as  the  “sailing  partials.” 
The  subscript  k  represents  the  leg  number  of  the  current  observation,  while  j  represents  the  leg 
number  of  a  particular  unknown.  In  each  evaluation  of  equation  (4),  k  remains  fixed  but  j  runs  over 
all  legs  in  the  problem. 

There  are  three  fundamental  points  to  be  noted  in  evaluating  equation  (4) : 

(1)  Each  partial  derivative  that  appears  in  equation  (4)  is  numerically  evaluated  for  the  time 
and  estimated  position  of  the  observation  being  processed,  using  the  current  best  estimates  of 
the  values  of  the  parameters  of  the  problem. 

(2)  When  the  subscripts  in  the  “numerator”  and  “denominator”  of  the  sailing  partials  are 
the  same,  that  is,  when  j  =  k,  the  partials  correspond  to  those  of  the  single-leg  case.  The 
expressions  from  the  appendix  to  [1]  can  be  used  directly,  although  the  notation  is  somewhat 
different.  For  example,  dg^jd^^  is  evaluated  using  the  equation  from  [1]  for  dg/d(po,  applied  to 
leg  2. 

(3)  The  summation  in  the  second  line  of  equation  (4)  effectively  ends  at  j  =  k,  since  all  the 
sailing  partials  are  zero  when  j  >  k,  for  reasons  stated  in  the  main  body  of  the  paper. 

The  major  issue  addressed  in  this  appendix  is  the  evaluation  of  the  sailing  partials  in  equation  (4) 
when  j  <  k.  These  occur  in  the  computation  of  the  coefficients  of  unknowns  relating  to  legs  preceding 
the  one  on  which  the  observation  was  taken.  To  facilitate  the  computation  of  these  partials,  we  will 
add  an  artificial  observation  at  the  end  of  each  leg;  that  is,  for  leg  j  the  added  computation  point  is 
at  t  =  tj+1.  We  compute  all  the  sailing  partials  that  appear  in  equation  (4)  for  this  point  just  as  if 
it  were  an  actual  observation  at  the  very  end  of  the  leg.  The  values  of  these  partials  are  stored  for 
future  reference.  This  point  is  not  included  in  the  solution  since  it  is  not  an  actual  observation  (we 
have  no  value  for  a^).  The  artificial  observation  is  processed  this  way  for  each  leg  except  the  final 
one,  where  it  is  superfluous.  The  utility  of  these  artificial  observations  will  soon  become  evident. 


First  Leg.  The  sailing  formulas  for  the  first  leg,  which  apply  to  times  t  where  t\  <  t  <  t2l  can 
be  represented  as 


Ht)  =  h  = 

^(t)  =  9i  =  9  (t  —  1 1 ,  (/>r ,  Ai ,  Ci ,  S\ ) 


(Al) 
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For  observations  on  the  first  leg,  where  k  =  1,  we  only  have  to  deal  with  j  =  1  since  the 
sailing  partials  are  all  zero  for  j  >  k.  Since  the  only  subscript  we  encounter  in  either  “numerators” 
or  “denominators”  is  1,  all  the  observations  on  the  first  leg  are  processed  just  as  in  the  single-leg 
problem. 

Second  Leg.  The  sailing  formulas  for  the  second  leg,  which  apply  to  times  t  where  t2  <  t  <  f3, 

are: 


4>{t)  —  fi  —  f  (t  -  t2l  <f>2,  a2,  c2l  s2) 

(A2) 

X(t)  —9 2  —  9  (t  —  1 2 ,  4>2 ,  A2  ,  C2  ,  S2  ) 

However,  there  is  a  constraint  on  the  point  ((p2,  A2),  which  is  the  start  of  leg  2:  it  is  also  the  end  of 
leg  1.  That  is,  at  t  =  t2,  we  have  /2  =  fy  and  g2  =  gy.  More  explicitly, 


(/)2  —  </>(t2)  —  /2  —  fl  —  f  (h  —  t\i  ^1)  Sy 

t2  1 2 


A2  —  A(t2)  —  92 


=  g  i 


—  9  (t2  —  1 1 ,  (/>i ,  Ai ,  Ci ,  Si 


(A3) 


where  the  vertical  bar  means  “evaluated  at.” 

For  observations  on  the  second  leg,  where  k  =  2,  we  have  to  deal  with  j  =  1  and  2;  the  sailing 
partials  are  all  zero  for  j  >  k. 

Partials  for  Afy  and  AXy.  To  obtain  expressions  for  the  sailing  partials  needed  in  the  first  line 
of  equation  (4),  we  substitute  the  expressions  for  <f>2  and  A2  from  the  right  side  of  equations  (A3) 
into  equations  (A2),  then  differentiate  using  the  generalized  chain  rule.  The  results  are: 


df2  _  djfdjf  df2  dgy 

8(f)y  df>2  d(f)y  d\2  d(f)y 

df2  _  df^dfy^  df2  dgy 

d\y  df>2  d\y  d\2  d\y 

(A4) 

%2  dfhdjx  %2  dgy 

d<t>y  d<t>2  d(f)y  d\2  d(f)y 

dg2  _  dfhdjf  dg2  dgy 

d\y  df>2  d\y  d\2  8Xy 

On  the  right  sides  of  equations  (A4),  the  partials  of  f2  and  g2  (df2/d(j)2,  etc.,  the  first  factors  in  each 
term)  are  evaluated  for  the  particular  time  of  the  observation  being  processed,  and  are  computed 
from  the  single-leg  formulas  just  as  if  leg  2  were  the  only  leg  in  the  problem.  However,  the  partials  of 
fy  and  gy  (dfy/dfy,  etc.,  the  second  factors  in  each  term)  are  evaluated  for  time  t  =  t2  (as  required 
by  equation  (A3)),  so  their  values  have  already  been  computed — they  are  the  partials  computed  for 
the  artificial  observation  at  the  end  of  leg  1. 

Partials  for  ACy  and  A Sy .  Similarly,  we  differentiate  equations  (A2) ,  with  <f>2  and  A2  substituted 
from  equations  (A3) ,  to  obtain  the  sailing  partials  needed  for  the  second  line  of  equation  (4)  when 
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j  =  1: 


<9/2 

_  dh_ 

dh 

+ 

df2 

dgi 

dC\  ' 

~  df2 

dC\ 

dX2 

dC\ 

d/2 

_  dh_ 

dh 

+ 

dh 

dgi 

dSi  ' 

~  df2 

dS1 

dX2 

dSi 

%2 

_  dg2 

dh 

+ 

dg2 

dgi 

dC\  ' 

~  df2 

dC\ 

dX2 

dC\ 

%2 

_  dfh 

dh 

+ 

dg2 

dgi 

dSi  ' 

~  d<j)2 

dS1 

dX2 

dSi 

(A5) 


The  partials  of  f2  and  g2  on  the  right  side  of  equation  (A5)  are  the  same  as  those  on  the  right  side 
of  equation  (A4).  The  partials  of  h  and  g1  are  evaluated  for  time  t  =  t2,  the  end  point  of  leg  1,  so 
their  values  have  also  been  previously  computed. 


Partials  for  AC2  and  A S2.  For  the  second  line  of  equation  (4)  when  j  =  2,  we  encounter  only 
partials  where  the  “numerator”  and  “denominator”  have  the  same  subscript  (i.e. ,  j  =  k  =  2).  These 
are  evaluated  for  each  observation  from  leg  2  just  as  if  it  were  a  single-leg  case. 


Third  Leg.  The  sailing  formulas  in  latitude  and  longitude  for  leg  3  can  be  represented  as 


4>(t)  —  /3  —  f  (t  -  t3l  <f> 3,  X3lC3l  S3) 
X(t)  =  93  =  9  (t  -  t3l  (j> 3,  X3lC3l  S3) 


(A6) 


We  have  the  usual  constraint  on  the  point  (c/>3,  A3),  which  is  the  start  of  leg  3,  but  also  the  end  of 
leg  2: 


4>  3  —  (t){t'i)  —  /3  ~  f2 

f3 


A3  —  A  (t3)  —  g3\ 


=  9  2 


For  observations  on  the  third  leg,  where  k  =  3, 


/  (t3  —  1 2 , 4>2 ,  a2  ,  c2 ,  s2 

9  (f3  —  f  2  >  ^2  >  A2  ,  C2 ,  ^2 
we  have  to  deal  with  j 


1,  2,  and  3. 


(A7) 


Partials  for  Afi  and  AA^  The  sailing  partials  needed  for  the  first  hne  of  equation  (4)  are 
obtained  by  differentiating  equations  (A6),  with  and  A3  substituted  from  equations  (A7).  With 
the  help  of  the  chain  rule,  the  following  results  are  obtained: 


df3  _  djfdjf  df3  dg2 

dfi  dfo  dfi  <9A3  dfi 

df3  _  djfdJji  df3  dg2 

d\\  dfo  d\i  <9A3  9Ax 

(A8) 

%3  dfpdh_  %3  dg2 

d(f)1  d(j>3  d(f)1  d\3  d(f)1 

dg3  _  dffcdh.  dg3  dg2 

d\\  d(f)3  d\i  dX3  dX1 

Flere,  the  partials  of  f3  and  g3  on  the  right  side  are  evaluated  for  the  particular  time  of  the  observation 
being  processed,  from  the  single-leg  formulas  applied  to  leg  3.  The  partials  of  f2  and  g2  are  evaluated 
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for  time  t  =  t3  ( a  requirement  of  equation  (A7)),  so  we  already  have  their  values — they  were 
computed  for  the  artificial  observation  at  the  end  of  leg  2  when  equation  (A4)  was  applied. 

Partials  for  A C\  and  A S\ .  The  sailing  partials  needed  for  the  second  line  of  equation  (4)  when 
j  =  1  are  obtained  by  the  same  procedure  as  used  for  equations  (A8): 


<9/3 

_  djf 

<9/2 

+ 

d/s 

dg2 

dC\  ~ 

~  df3 

dC\ 

d\3 

dC\ 

d/3 

_  d/s 

<9/2 

+ 

d/s 

dg2 

dSi  ~ 

~  df3 

dSi 

d\3 

dSi 

%3 

_  %3 

<9/2 

+ 

dg3 

dg2 

dC\  ~ 

~  df3 

dC\ 

d\3 

dC\ 

%3 

_  %3 

<9/2 

+ 

dg3 

dg2 

dSi  ~ 

~  df3 

dSi 

d\3 

dSi 

(A9) 


The  partials  of  f3  and  g3  on  the  right  are  the  same  as  those  in  equation  (A8).  The  partials  of  f2 
and  g2  are  evaluated  for  time  t  =  t3,  the  end  point  of  leg  2,  and  their  values  were  computed  when 
equation  (A5)  was  evaluated  for  the  artificial  observation  at  the  end  of  leg  2. 

Partials  for  AC2  and  A S2.  The  sailing  partials  used  in  the  second  line  of  equation  (4)  when 
j  =  2  are: 
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+ 
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dg2 

dC2  ' 

~  df3 
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dfa 

dC2 
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dg3 

dg2 
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~  df3 

dS2 

d\3 

dS2 

The  partials  of  /3  and  g3  on  the  right  side  of  equation  (A10)  are  the  same  as  those  in  equations  (A8) 
and  (A9).  The  partials  of  f2  and  g2  are  evaluated  for  time  t  =  t3,  the  end  point  of  leg  2,  and  their 
values  are  also  available  from  previous  computation. 

Partials  for  A  C3  and  A  S3.  For  the  second  line  of  equation  (4)  when  j  =  3,  all  the  sailing 
partials  have  the  same  subscript,  3,  in  their  “numerator”  and  ’’denominator.”  These  are  evaluated 
for  the  observation  using  the  single-leg  formulas  applied  to  leg  3. 

General  Strategy.  The  pattern  should  be  obvious  by  now.  The  partials  of  the  sailing  formulas 
that  appear  in  equation  (4)  have  the  general  form  dy k  / dxj ,  where  y  represents  either  the  function 
/  or  <7;  x  represents  any  of  the  four  parameters  <f>,  A,  C,  or  S';  k  is  the  number  of  the  current 
observation’s  leg;  and  j  is  an  index  that  represents  any  of  the  N  legs  in  the  problem.  When  j  >  k 
the  partial  is  zero.  The  non-zero  partials  are  of  two  types.  In  the  first  type,  j  =  k,  that  is,  the 
partial  has  the  form  dyk/dxi~.  These  partials  are  computed  from  the  single-leg  formulas  given  in 
[1],  applied  to  leg  k,  and  are  evaluated  for  the  time  of  the  observation  being  processed.  The  second 
type  of  partial  has  j  <  k,  and  we  have  seen  that  these  partials  are  expanded  as  follows: 


dyk  _  dyk  dfk- 1  dyk 

dxj  d(f)k  dxj  d\k  dxj 


(All) 
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Looking  at  the  two  terms  on  right  side  of  equation  (All),  we  see  that  the  first  factor  in  each  is  a 
partial  of  the  first  type,  which  we  know  how  to  evaluate.  For  a  given  observation,  these  first  factors 
are  used  repeatedly.  The  second  factor  is  a  partial  that  must  be  evaluated  for  time  t  =  t^,  which  is 
the  instant  at  the  end  of  leg  k  —  1  (the  previous  leg) .  Since  we  have  inserted  an  artificial  observation 
at  the  end  of  each  leg,  these  partials  have  already  been  computed.  For  each  sailing-formula  partial 
that  appears  in  equation  (4),  expanded  as  above,  these  second  factors  remain  the  same  for  all 
observations  within  a  leg.  It  should  be  noted  that  some  of  the  partials  used  in  equation  (All)  are 
zero  if  the  rhumb-line  sailing  formulas  from  [1]  are  used.  (For  example,  latitude  is  not  a  function  of 
longitude,  so  dfk/d A*,  =  0.)  The  algorithm  has  been  presented  in  its  fullest  generality  so  that  it  is 
not  limited  to  a  specific  set  of  sailing  formulas. 

This  algorithm  can  be  considered  an  extension  of  the  one  used  for  the  single-leg  case.  It  requires 
no  new  mathematics  beyond  trivial  multiplications  and  additions.  Equation  (All)  amounts  to  a 
recurrence  relation  that  builds  the  sailing  partials  for  the  j  <  k  case  from  those  already  computed  for 
the  previous  leg,  together  with  partials  computed  using  the  single-leg  formulas  applied  to  the  current 
leg.  No  changes  in  the  least-squares  procedure  are  needed  beyond  providing  for  the  extra  unknown 
parameters  that  the  multiple-leg  case  requires.  Using  this  algorithm,  a  computer  implementation  of 
the  single-leg  case  can  be  modified  to  handle  the  multiple-leg  case  with  relatively  little  difficulty — 
the  additional  code  is  mostly  involved  with  bookkeeping.  Some  new  arrays  are  required  for  storage 
of  the  sailing  partials  computed  for  the  artificial  observations  at  the  end  of  each  leg.  However,  the 
storage  required  is  small  since  in  practice  the  number  of  legs  that  would  be  included  in  a  solution 
would  probably  never  exceed  about  five. 

As  mentioned  in  the  main  text,  for  practical  applications,  the  algorithm  should  in  most  cases  be 
implemented  in  a  time-reversed  sense.  That  is,  leg  1  should  be  the  most  recent  (latest)  leg  and  leg  N 
should  be  the  first  (earliest)  leg.  That  allows  the  fix  (</>i,Ai),  which  can  be  anywhere  on  leg  1,  to 
represent  the  current  position  of  the  observer’s  vessel  if  desired.  In  a  time-reversed  implementation, 
the  observations  are  taken  in  reverse  chronological  order.  Each  of  the  artificial  observations  must 
then  be  added  at  the  beginning  of  its  leg,  so  that  it  is  the  last  computational  point  encountered  on 
the  leg  as  the  observations  are  processed. 
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Appendix  B 

Alternative  Multiple-Leg  Solution  Algorithms 


What  alternative  algorithms  are  available  for  the  multiple-leg  navigational  solution?  As  men¬ 
tioned  in  the  main  part  of  the  paper,  the  multiple-leg  case  is  an  example  of  a  constrained  least-squares 
problem.  Least-squares  theory  allows  constraint  equations — known  relationships  among  the  param¬ 
eters  being  solved  for — to  be  included  as  part  of  the  solution  (see,  for  example,  [8],  Chapter  9). 
In  principle,  this  theory  would  allow  us  to  use  equation  (3),  with  four  unknown  parameters  per 
leg,  in  combination  with  linearized  versions  of  constraint  equations  such  as  (A3)  and  (A7)  in  Ap¬ 
pendix  A.  However,  this  approach  is  considerably  more  difficult  computationally,  involving,  among 
other  things,  another  matrix  inversion.  Generally,  “canned”  least-squares  packages  do  not  provide 
for  this  type  of  solution  (although  at  least  one,  IMSL  [10],  does). 

Another  approach,  discussed  in  [8],  uses  the  constraint  equations  to  eliminate  the  dependent 
(redundant)  parameters  in  the  problem.  That  is,  the  constraints  are  used  to  solve  for  a  set  of 
independent  (free)  parameters  “up  front,”  and  all  the  relevant  equations  are  rewritten  in  terms 
of  this  set  of  independent  unknowns.  The  algorithm  outlined  above  is  actually  a  variant  of  this 
approach.  In  our  case,  the  coordinates  of  the  starting  points  of  legs  2  through  N  are  the  dependent 
parameters.  The  algorithm  does  not  explicitly  solve  for  and  eliminate  these  parameters — rather,  the 
algorithm  works  entirely  with  the  derivatives  of  the  parameters,  and  uses  the  chain  rule  to  effectively 
remove  the  starting  point  coordinates  of  legs  2  through  N  from  the  solution.  The  solution,  therefore, 
involves  only  a  set  of  independent  parameters. 

The  direct  elimination  of  the  dependent  parameters  is  difficult  with  the  sailing  formulas  used 
above  but  the  method  can  be  illustrated  if  the  vessel’s  track  is  approximated  by  piecewise  polyno¬ 
mials.  (Spline  theory  is  not  really  applicable,  since  continuity  in  derivatives  is  not  required  at  the 
junction  points.)  For  short  rhumb  lines,  latitude  can  be  represented  as  a  linear  function  of  time, 
although  longitude  requires  a  quadratic  polynomial  (great  circles  would  require  quadratics  in  both 
coordinates).  Suppose  we  look  at  the  equations  for  longitude  for  a  three-leg  problem: 

Leg  1:  A  (f)  =  g\  =  Ax  +  Ax  (t  —  H)  +  \^\{t  —  G)2 

Leg  2:  \(t)  =  g2  =  \2  +  \2(t  —  t2)  +  A2(f  —  t2)2  (Bl) 

Leg  3:  \(t)  =  g3  =  A3  +  \3(t  —  t3)  +  ^  A3(f  —  t3)2 

These  equations  contain  nine  parameters:  Ai,  Ai,  Ai,  A2,  A2,  A2,  A3,  A3,  and  A3.  However,  two  of 
these,  A2  and  A3,  the  longitudes  at  the  beginning  of  legs  2  and  3,  are  dependent  parameters.  They 
can  be  related  to  the  other  parameters  by  the  usual  constraints,  which  provide  for  continuity  in 
position  at  the  leg  intersections.  For  this  case,  the  constraints,  analogous  to  equations  (A3)  and 
(A7)  in  Appendix  A,  are  g2  =  g1  at  t  =  t2  and  g3  =  g2  at  t  =  t3.  These  conditions  allow  us  to 
eliminate  A2  and  A3  from  equations  (Bl),  yielding 


91  =  Ai  +  Ai  (t  —  t1)  +  ^A1(t  —  t])2 

92  =  +  ^i(G  —  H)  +  5^1  (G  —  H)2  +  ^2(f  —  G)  +  5^2(f  —  G)2 

93  =  +  ^i(G  —  H)  +  (G  —  H)2  +  ^2(^3  —  t2)  +  ^2(^3  —  G)2 

+  A3(f  —  t3)  +  ^\3(t  —  t3)2 


(B2) 
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Equations  (B2)  are  the  new  sailing  formulas  for  longitude  to  be  used  for  the  three  legs,  which 
involve  no  dependent  parameters.  The  latitude  equations  (which  have  no  quadratic  terms)  can 
be  treated  similarly.  We  are  left  with  only  one  geographic  point  in  the  problem,  defined  by  the 
parameters  (pi  and  Ai,  plus,  for  each  leg  j,  three  motion  parameters,  <f>j,  A j,  and  A j.  We  can 
then  form  a  conditional  equation,  similar  to  equation  (4)  in  the  main  text,  but  with  these  eleven 
parameters  as  unknowns.  A  least-squares  solution  would  then  be  completely  straightforward. 

This  seems  to  be  a  simpler  approach  to  the  problem,  but  there  are  prices  to  be  paid.  First, 
the  motion  parameters  for  each  leg  j — that  is,  <f>j,  A  j,  and  A  j — are  somewhat  abstract  and  would 
have  to  be  transformed  to  quantities  such  as  course  or  speed  that  would  be  useful  to  the  navigator. 
Second,  there  are  now  three  motion  parameters  to  solve  for  per  leg  instead  of  two,  requiring  more 
observations  for  a  reliable  solution.  Third,  the  polynomial  representation  of  each  leg’s  track  is 
simply  a  Taylor  expansion  of  more  exact  formulas  (here,  for  rhumb  lines)  and  extra  terms  might 
be  required  for  long  legs,  depending  on  the  accuracy  requirements.  Every  extra  term  adds  another 
free  parameter  to  the  problem,  and  the  requirement  for  more  observations.  The  algorithm  presented 
in  the  main  text  of  this  paper  and  Appendix  A  is  thus  more  economical  in  its  use  of  the  available 
observations,  a  critical  consideration. 
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